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Abstract 

We study exchange of stability in the dynamics of solitary wave 
solutions under changes in the nonlinear balance in a 1+1 evolu- 
tionary partial differential equation related both to shallow water 
waves and to turbulence. We find that solutions of the equation 
m t + um x + bu x m — vm xx with m = u — a 2 u xx for fluid velocity 
u(x, t) change their behavior at the special values 6 = 0, ±1, ±2, ±3. 

PACS numbers: ll.10.Lm, 03.50.-z, 05.45.Yv, 05.45.-a 
Keywords: solitary waves, solitons, peakons, nonlinear evolution, tur- 
bulence modeling 
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1 Introduction 

1.1 Linear and nonlinear balances in shallow water waves 

The primary physical mechanism for the propagation of solitary shallow wa- 
ter waves is the balance between nonlinear steepening and linear dispersion. 
This balance appears uniquely at linear order in an asymptotic expansion in 
the Korteweg-de Vries equation (KdV), 

3ei 3e2 

u t + c u x + —uu x + — u xxx = . (1) 

Here the expansion parameters satisfy e\ > €2 > ef and are defined by e± = 
a/h and €2 = h 2 /l 2 , in terms of a, h and I, which denote wave amplitude, 
mean water depth, and typical horizontal length scale (e.g., a wavelength), 
respectively. KdV possesses the famous sech 2 solitary wave solution u(x, t) = 
uq sech 2 ((x — ct)y/uo/j/2), for no = 2(c — co)/ei and 7 = 3e2/(5ei), see pQ. 
The Benjamin-Bona-Mahoney equation (BBM), 

3ej 3e2 1 , . 

ik + c u x + — uu x - — c u xxt = , (2) 

has the same sech 2 shape for its solitary wave, but with 7 replaced by 
7' = C7/C0. BBM is asymptotically equivalent to KdV at order 0(61,62). 
However, the linear dispersion relation for BBM is a better match than that 
for KdV in the physical application to shallow water waves. 

Beyond KdV at linear order, the asymptotic expansion at quadratic 
order in the small parameters e\ and £2 does not produce a unique wave 
equation for unidirectional shallow water waves |2j. Instead, the asymptotic 
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expansion produces an entire family of shallow water wave equations that 
are asymptotically equivalent to each other at quadratic order in the shal- 
low water expansion parameters. The equations in this family are related 
amongst themselves by a continuous, three-parameter group of nonlinear, 
nonlocal transformations of variables given by, 



u 



v + ei(aiv 2 + a 2 v x d x 1 v) + e 2 a 3 v xx , (3) 



in which (01,02,03) are real parameters. This transformation group was 
first introduced for determining normal forms of shallow water equations by 
Kodama in pj. 

Among the family of asymptotically equivalent shallow water wave equa- 
tions at quadratic order accuracy in the small parameters e\ = a/h and 
€2 = h 2 /l 2 are several equations that are completely integrable. As for KdV 
at linear order, these integrable shallow water equations at quadratic order 
possess soliton solutions that interact via elastic collisions. In particular, the 
equation in the KdV hierarchy with fifth-order derivatives is present among 
these integrable equations, as shown in [I]. 

The family of shallow water equations at quadratic order accuracy that 
are asymptotically equivalent under Kodama transformations (jHJ) contains 
the following sub- family, derived in [2J, 

362 

m t + c u x + —u xxx + ei(um x + bmu x ) = 0, (4) 
Evolution Digpersion Nonlinearity 

where m = u — (19e2/G0)u xx . The set of 1+1 evolutionary equations ((3J) is 
nonlocal, dispersive and nonlinear. The nonlinearity in this set of asymptot- 
ically equivalent shallow water equations is parameterized by the real con- 
stant 6, which depends on the group parameters (ai, 02, 03) in the Kodama 
transformation ©. An asymptotically equivalent shallow water equation 
for any 6 7^ — 1 may be achieved by a Kodama transformation. However, 
the case b = — 1 violates the asymptotic ordering and the corresponding 
Kodama transformation is singular in this case [2j. 

The cases 6 = 2 and 6 = 3 are special values for the 6— equation The 
case 6 = 2 restricts (J3} to the integrable Camassa-Holm equation (CH) 0. 
The case 6 = 3 in @ is the Degasperis-Procesi equation (DP) [B], which 
was shown to be integrable in jj]. These two cases exhaust the integrable 
candidates for Q , as was shown using Painleve analysis in [7j . The b-family 
of equations (|1J) was also shown in [H] to admit the symmetry conditions 
necessary for integrability only in the cases 6 = 2 for CH and 6 = 3 for DP. 
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KdV (UJ) and the cases 6 = 2 (CH) and 6 = 3 (DP) of equation (gj) are 
three completely integrable Hamiltonian equations that possess solitons as 
traveling waves. In KdV, the balance that confines the traveling wave soli- 
ton occurs between nonlinear steepening and linear dispersion. This is the 
leading order asymptotic balance for shallow water waves. However, even 
in the absence of linear dispersion, the parameter 6 in equation © intro- 
duces additional possibilities for balance, including the nonlinear /nonlocal 
balance in the following (rescaled) dispersionless case of CH that was studied 
previously for 6 = 2 in |Hj , 

mt + um x + 6 mu x = , with m = u — o?u xx . (5) 

The nonlinear /nonlocal balance in this equation, even in the absence of 
linear dispersion, can still produce a confined solitary traveling wave pulse 
u(x,t) = ce~\ x ~ ct \l a , called the peakon [5]. The peakon moves with speed 
equal to its amplitude and has a jump in derivative at its peak. Peakons 
for either 6 = 2 or 6 = 3 are true solitons that interact via elastic collisions 
under CH dynamics, or DP dynamics, respectively. In addition, the CH 
and DP initial value problems are both completely integrable as Hamilto- 
nian systems by using the inverse spectral transform (1ST) method for an 
isospectral linear eigenvalue problem whose purely discrete spectrum gives 
the asymptotic speeds of the peakons [7]. Figure ^ shows the evolution 
under dispersionless CH for the case 6 = 2 in equation (0 ) of a Gaussian 
initial velocity distribution of unit area and width 5a. Peakon solutions 
exist for equation @ with any value of 6. However, we shall find that the 
stability of these peakon solutions requires 6 > 1. 

The properties of dispersionless CH for the case 6 = 2 in equation © 
and also the related class of dispersionless equations consisting of © with 
u = g * m = g(x — y)m(y)dy for the case 6 = 2 and an even kernel 
g(x) = g(—x) including g{x) = e~\ x ^ a are studied further in Fringer and 
Holm 9 . The properties of dispersionless DP for the case 6 = 3 in equation 
(0) are studied further in [7j, [TUj . 

The dispersionless limit of KdV in upon rescaling velocity u, yields 
the Burgers equation 

u t + uu x - v u xx = , (6) 

in which we now also add constant viscosity v. The characteristic Burgers 
solution is the classic ramp and cliff shown arising from the Gaussian initial 
condition in Figure 121 In the ramp/cliff solution, nonlinear steepening is 
balanced by linear viscosity to produce the "cliff" whose width is controlled 
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Figure 1: Evolution under equation © with b = 2 of a Gaussian initial 
velocity distribution of unit area and width 5a. 

by the magnitude of viscosity v. The "ramp" is the self-similar u ~ x/t part 
of the solution for which the viscous term vanishes. 

The equations KdV, BBM, CH, DP, the other b— equations and Burgers 
all preserve the area M = J^*L u dx (total momentum) for a solution that 
vanishes at spatial infinity. 

1.2 Problem statement for the present work 

In this paper, we shall treat the constant b in the dispersionless case of 
equation Q) as a bifurcation parameter. We shall study numerically the 
exchange of stability among three types of solitary wave solutions, as b is 
varied. This is done by using numerical simulations to investigate its initial 
value problem on a periodic domain of length L»a for various values of 
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Figure 2: The Burgers solution is the classic ramp/cliff shown arising from 
the Gaussian initial condition. 

the parameter 6. When necessary, we shall also add viscosity in the form 1 

mt + um x + bu x m = vm xx , (7) 
convection stretching viscosity 

with m = u — cP"u xx . We shall treat the following cases, which we shall show 
analytically are special values: 6 = 0, dbl, ±2, ±3. 



2 Exchange of stability 

2.1 Numerics for exchange of stability from b = to b = ±2 

We begin by considering the initial value problem for case b = in equation 
(|7)l. The first part of Figure 01 shows for 6 = that a peakon of width 
w = 5a is unstable to forming a ramp/cliff. Namely, the initial peakon 
develops into a Burgers-type ramp/cliff solution for 6 = 0, just as it would 

1 Relation of the 6— equation (J4j| to 3D turbulence closure models. Higher- 
dimensional representatives of the viscous b— family of equations J3J have appeared re- 
cently in 3D turbulence modeling. The b = case of equation Q is the ID version of 
Leray's regularization of the Navier-Stokes equations in 3D Leray regularization has 
recently been revived as an approach to deriving turbulence closure models for large eddy 
simulations of turbulent mixing layers in 3D For 3D studies of the turbulence model 
corresponding to the b~ equation J1J with 6 = 2, see |13|. |14|. |15|. 
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Figure 3: Exchange of stability between the classic ramp/cliff for 6 = and 
the peakons for 6 = 2. 

have done for the Burgers equation. Note, however, that the case 6 = in 
equation is not the Burgers equation. At time T = 2500, we stopped this 
part of the numerical calculation and used the result as the initial condition 
for continuing the 6— equation evolution with 6 changed to 6 = 2, a case 
that supports peakons. The remainder of Figure shows the exchange of 
stability after 6 = is changed to 6 = 2. Thus, the Burgers-type ramp/cliff 
solution is stable for 6 = in equation Q. However, when treated as an 
initial condition, the ramp/cliff solution is unstable and breaks into a train 
of peakons for 6 = 2. Figure0]shows the spatial profiles of these solutions at 
the initial time, at the end of the 6 = evolution, and at the end of the 6 = 2 
evolution. A similar exchange of stability occurs between the ramp/cliff for 
6 = and the train of peakons for 6 = 3. 

If we switch from 6 = to 6 = — 2 instead of 6 = 2, then Figure El shows 
that the ramp/cliff at time T = 2500 evolves into a train of leftward moving 
solitary waves (leftons). Thus, a different exchange of stability occurs be- 
tween ramp/cliff solutions for 6 = and leftons for 6 = —2 (and for 6 = —3 
as well). Figure IHl plots the comparison between the stationary solution 
u = uq sech 2 (x/(2a)) for 6 = —2 and the final states for an initial Gaussian 
of width w = 10 evolved to leftons. Similar agreement occurs for the case 
6 = —3, whose stationary solution is u = uq sech (x/a). Thus, confined 
initial states move leftward and evolve into stationary solutions for 6 = — 2 
and 6 = —3. 
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Figure 4: Spatial profiles of the solutions at the initial time and at the end 
times of the 6 = evolution and the 6 = 2 evolution. 
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Figure 6: Comparison between leftons and the stationary solution u = 
uq sech 2 (x/(2a)) for 6 = —2. 

2.2 Analytical discussion of the exchange of stability 

The exchange of stability as the parameter 6 varies can be understood ana- 
lytically as arising from the following properties of the b— equation l|Tj). 

• For < b < 1, the L 1 / 6 norm of the momentum 1 1 rrx 1 1 = ^ f^°oo {ml^j 
is preserved by Q, in the inviscid case for which v = 0. Thus, the ramp/cliff 
solution for 6 = forms in a regime of classical analytical control in which 
\m\ possesses a maximum principle. Numerically, we observe: (1) the cliff's 
width is controlled by the magnitude of a; and (2) the viscosity v can be 
made negligible without measurably affecting the computed solution. 

• For 1 < b < 3 in the inviscid case, an inflection point with negative 
slope may be shown to develop into a vertical slope in finite time. However, 
numerically, we observe the nonlinear evolution conspires to eliminate such 
singular inflection points, e.g., by forming peakons, which are weak solutions. 

• For b = — 1 the lowest order nonlinearity in equation ((TJ) (the steepening 
term uu x ) has coefficient zero. The evolution then develops from the higher 
order nonlinear terms. This can be seen by rewriting the b— equation (J7J) 
equivalently as 



ut + (b + l)uu x — v u 



C? (lit + uu : 



'X 



6 — 3 2 \ 
vu xx + —^-u x ) 



(8) 



Thus, solutions with small curvature are nearly stationary for 6 = 



-1. The 
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leading order nonlinearity is absent, so the dynamics is governed by the 
balance between viscous diffusion and higher-order terms. The case a — > oo 
in (JHJ) is interesting, as well. However, this case will be discussed elsewhere. 
• Equation (J7| has stationary (c = 0) solutions u = uq sech 2 (x/(2a)) for 
b = — 2 and u = itosech(x/a) for b = —3. Figure El shows that the evolution 
of equation ((7J) rapidly approaches the stationary solution for b = —2. (A 
similar result holds for b = —3.) The linearized evolution of (JZJ) around these 
stationary solutions is stable for all wavelengths. Therefore, for any nonzero 
viscosity, these stationary solutions have a nontrivial basin of attraction. 
Note: equation (jJJ) with either b = —2, or b = —3 is not Hamiltonian. The 
mathematical properties of the b— equation in these cases is discussed in |10| . 

2.3 Confirmations of our numerical methods 

In our numerical runs we advanced the b— equation (jJJ) with an explicit, 
variable timestep fourth/fifth order Runge-Kutta-Fehlberg (RKF45) predic- 
tor/corrector. We selected the timestep for numerical stability by trial and 
error, while our code selected the timestep for numerical accuracy (not to 
exceed the timestep for numerical stability). We used a very strict relative 
error tolerance per timestep, e = 10 -8 . 

We computed spatial derivatives using 4th order finite differences, gen- 
erally at resolutions of 2 13 = 8192 or 2 14 = 16384 zones. To invert the 
Helmholtz operator in determining u(x, t) from m(x, t), we used the Fourier 
transform. When the numerical approximation of the nonlinear terms had 
aliasing errors in the high wavenumbers, we removed these errors by apply- 
ing a high pass filtered artificial viscosity. 

The quality of the numerical convergence may be checked analytically in 
the collision of two peakons of asymptotic speeds c\ and C2, for which the 
minimum peakon separation with b > 1 is 



For 6 = 2, c\ = 1, p2 = 1/2 and a = 5, formula © implies q m in = 
10 In 3 = 10.9861. Our numerical results with the resolution of 2 14 zones 
yield q m in = 11.0049. The small discrepancy, less than 0.2%, occurs largely 
because our numerical measurement of q m i n is obtained by examining the 
peakon positions at each internal timestep in the code without interpolation, 
although one is unlikely to land exactly on the time at which the minimum 
separation occurs. The code's true accuracy is better than the above mea- 
sure indicates, because the intermediate steps involved in advancing the 




(9) 



D. D. Holm <fe M. F. Staley Solitons, peakons, ramps/cliffs and leftons 11 



solution from one discrete time to the next with an RKF45 method cancel 
the higher-order discretization errors. 

Likewise, for peakons with 6 = 3, formula © gives the minimum separa- 
tion when ci = 1, C2 = 1/2, and a = 5, as q m i n = 5 In (3/ (3— VS)) = 14.3068. 
In this case, our numerical results yield q m in = 14.2924, a discrepancy of 
only 0.1%. 




X 



Figure 7: Scores of peakon collisions arise for 6 = 3 (DP) starting from an 
initial Gaussian velocity distribution of width w = 10. 

Of course, the two-peakon collision is rather simple compared to the 
plethora of other multi-wave dynamics that occurs in this problem, as in 
Figures [7| and |HJ For this reason, we also checked the convergence of our 
numerical algorithms by verifying that the relative phases of the peakons in 
the various figures remained invariant under grid refinement. Moreover, the 
integrity of the waveforms in our figures attests to the convergence of the nu- 
merical algorithm - after scores of collisions, the peakon waveforms are still 
extremely well preserved. The preservation of these peakon/soliton wave- 
forms after so many collisions would not have occurred unless the numerics 
had converged well. 

3 Conclusions and further comments 

Equation (J5J) comprises a family of reversible, parity invariant, evolutionary 
1+1 PDEs that arise in the zero-dispersion case of shallow water waves at 
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Figure 8: Under DP in © with 6 = 3, many peakons of unit width for 
a = 1 emerge from an initial pair of peakons of width w = 5. 

quadratic order in the standard asymptotic expansion, 

mt + um x + b mu x = , with m = u — a 2 u xx . (10) 

The paper identified the bifurcations and exchanges of stability of its trav- 
eling wave solutions as a function of the nonlinear balance parameter 6. 

Numerical computations of the initial value problem for (jlUj) when b > 1 
showed the emergence of its stable particle-like peakon solutions and their 
interactions. These may be obtained analytically by superposing N peakon 
traveling wave solutions for the inviscid case, u(x,t) = ce -\ x ~ ct \ as 

N N 

u(x,t) = Y,Pi(t)e~ l ' X ~ 9l(t)l and m(x, t) = ^ Pi {t)8{x - , (11) 
i=l i=l 

for any real constant b. For any b > 1, the peakons are stable and undergo 
particle- like dynamics in terms of the moduli variables Pi(t) and with 
i = 1,...,N. The peakon dynamics studied for b > 1 in this framework 
displayed all of the classical soliton interaction behavior for peakons found 
in H2 for the case 6 = 2. This behavior included pairwise elastic scat- 
tering of peakons, dominance of the initial value problem by confined pulses 
and asymptotic sorting according to height - all without requiring complete 
integrability. Thus, the "emergent pattern" for 6 > 1 in the nonlinear evo- 
lution governed by the 6— equation H10|) was the rightward moving peakon 
train, ordered by height. 
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A second type of emergent pattern of the initial value problem for (|1U[) 
occurs in the parameter region < b < 1. This is the classic Burgers- 
type ramp/cliff structure as in Figure |1J In contrast, for the parameter 
region b < — 1 a third type of behavior arises, consisting of leftward moving 
structures as in Figure El That the solution behavior should depend on the 
value of b is clear from the velocity form of equation ()10|) written in (|SJ) . 

Three regions of b. We found that the solution behavior for the 
b— equation (|1U|) changes its character near the two boundaries, b = ±1, 
of the following three regions in the balance parameter b. 

(Bl) In the stable peakon region b > 1, the Steepening Lemma for the 
b— equation (|10|) proven in [2] for 6 = 2 and in |10| for 1 < b < 
3 allows inflection points with negative slopes to escape verticality 
by producing a jump in spatial derivative at the peak of a traveling 
wave that eliminates the inflection points altogether. Peakon behavior 
dominates this region. When b < 1 we found the solution behavior of 
the b— equation changed its character and excluded the peakons 
entirely. 

(B2) In the Burgers region < b < 1, the L 1 / 6 norm of the variable m is 
controlled 2 and the solution behavior is dominated by ramps and cliffs, 
as for the usual Burgers equation. Similar Burgers- type ramp/cliff 
solution properties hold for the region — 1 < b < 0, for which the 
norm of the variable l/|m| is controlled. At the boundary of the 
latter region, for b = — 1, the active transport equation (jlOj) admits 
stationary plane waves as exact nonlinear solutions. However, shallow 
water asymptotic ordering is broken at b = — 1. 

(B3) In the steady pulse region b < — 1, pulse trains form that move left- 
ward from a positive velocity initial condition (instead of moving right- 
ward, as for b > —1). These leftward-moving pulse trains are found to 
approach a steady state. 

At the two boundaries, b = ±1, of these three regions in b, the emergent 
solutions of the b— family of equations are observed numerically to ex- 
change stability. Although we discussed some analytical indications, a full 
explanation of these properties remains open for future research. 

2 For b = 0, this is a maximum principle for \m\. 



D. D. Holm <fe M. F. Staley Solitons, peakons, ramps/cliffs and leftons 14 



Acknowledgments 

We are grateful to A. Degasperis, A. N. W. Hone, J. M. Hyman, S. Kurien, 
C. D. Levermore, R. Lowrie and E. S. Titi for their thoughtful insights and 
remarks. This work was supported by the US DOE under contracts W-7405- 
ENG-36 and the Applied Mathematical Sciences Program KC-07-01-01. 

References 

[1] M.J. Ablowitz and P.A. Clarkson, Solitons, Nonlinear Evolution Equa- 
tions and Inverse Scattering, Cambridge University Press (1991). 

[2] H.R. Dullin, D.D. Holm and G. Gottwald, Phys. Rev. Lett, 87, 194501 

(2001) . 

[3] Y. Kodama, Phys. Lett. A 107, 245-249 (1985); Ibid 112, 193-196 
(1985); Ibid 123, 276-282 (1987). 

[4] Z. Li and N.R. Sibgatullin, J. Appl. Maths. Mechs. 61, 177-182 (1997). 

[5] R. Camassa and D.D. Holm, Phys. Rev. Lett. 71 1661 (1993). 

[6] A. Degasperis and M. Procesi, Asymptotic integrability, in Symmetry 
and Perturbation Theory, edited by A. Degasperis and G. Gaeta, World 
Scientific (1999) pp. 23-37. 

[7] A. Degasperis, D. D. Holm and A. N. W. Hone, A new integrable equa- 
tion with peakon solutions. Theoret. and Math. Phys. 133, 1463-1474 

(2002) . |http://xxx.lanl.gov/abs/nlin.SI/0205023| 

[8] A.S. Mikhailov and V. Novikov, (private communication). 

[9] O. Fringer and D.D. Holm, Physica D 150 237 (2001). 

[10] D.D. Holm and M.F. Staley, Wave structures and nonlinear balances in 
a family of 1+1 evolutionary PDEs. To appear in SIAM J. Dyn. Syst. 

(2003) . 

[11] J. Leray, Sur le mouvement d'un liquide visqueux emplissant l'espace, 
Acta Math. 63, 193-248 (1934). Reviewed, e.g., in P. Constantin, C. 
Foias, B. Nicolaenko & R. Temam, Integral manifolds and inertial mani- 
folds for dissipative partial differential equations. Applied Mathematical 
Sciences, 70, (Springer- Verlag, New York-Berlin, 1989). 



D. D. Holm <fe M. F. Staley Solitons, peakons, ramps/cliffs and leftons 15 



[12] B.G. Geurts and D.D. Holm, Alpha-modeling strategy for LES of tur- 
bulent mixing. In Turbulent Flow Computation, edited by D. Drikakis 
and B. G. Geurts, Kluwer: London, pp. 237-278. 

[13] S.Y. Chen, C. Foias, D.D. Holm, E.J. Olson, E.S. Titi and S. Wynne, 
Phys. Rev. Lett. 81 5338 (1998). Ibid, Physica D 133 49 (1999). Ibid, 
Phys. Fluids 11 2343 (1999). 

[14] S.Y. Chen, D.D. Holm, L.G. Margolin and R. Zhang, Physica D 133 
66 (1999). 



[15] C. Foias, D.D. Holm and E.S. Titi, Physica D 152 505 (2001). Ibid, 
Dyn. and Diff. Eqns. 14 1 (2002). 



